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Abstract. We report on a number of careful numerical experiments motivated by the 
semiclassical (zero-dispersion, e J, 0) limit of the focusing nonlinear Schrodinger equation. 
Our experiments are designed to study the evolution of a particular family of perturbations 
of the initial data. These asymptotically small perturbations are precisely those that result 
from modifying the initial-data by using formal approximations to the spectrum of the 
associated spectral problem; such modified data has always been a standard part of the 
analysis of zero-dispersion limits of integrable systems. However, in the context of the 
focusing nonlinear Schrodinger equation, the ellipticity of the Whitham equations casts 
some doubt on the validity of this procedure. To carry out our experiments, we introduce 
an implicit finite difference scheme for the partial differential equation, and we validate both 
the proposed scheme and the standard split-step scheme against a numerical implementation 
of the inverse scattering transform for a special case in which the scattering data is known 
exactly. As part of this validation, we also investigate the use of the Krasny filter which is 
sometimes suggested as appropriate for nearly ill-posed problems such as we consider here. 
Our experiments show that that the 0{e) rate of convergence of the modified data to the 
true data is propagated to positive times including times after wave breaking. 



1. Introduction 

1.1. Background. In their pioneering work. Lax & Levermore [30 l [3T 1 132] used the inverse 
scattering transform (1ST) to study the zero-dispersion hmit of the initial- value problem for 
the Korteweg-de Vries equation: 

dty - Gyd^y = e^d^y , (1.1a) 

y{x,0) = yo{x). (1.1b) 

That is, they were able to characterize the limiting behavior of the family, indexed by e > 0, 
of solutions y^'^\x,t) of fll.lj) in the limit e J, 0. Their results (and those of others who have 
since extended and refined the analysis of fll.lj) . e.g., |^ W2\ \T5 \ \T2\ ) show — as one might 
guess — that for small times, y^^\x,t) converges strongly to y{x,t), the solution of 

dty-6ydxy = 0, (1.2a) 
y{x,0) = yo{x). (1.2b) 

For any fixed e > 0, solutions of (11. ip with smooth, decaying data exist and remain smooth 
for all if: > 0. By contrast, the limiting equation fll.2al) is known to support solutions which 
develop shocks in finite time regardless of the smoothness of the data. Indeed, a major 
impetus for the study of this problem has been to understand how the dispersive term in 
(11. lap prevents shock formation. Roughly, the solution develops rapid nonlinear oscillations 
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which carry energy away from a developing shock, and Lax-Levermore theory provides a 
rather precise description of the character of these oscillations in various regions of the xt- 
plane. An integral component of this description are the Whitham or modulation equations; 
these partial differential equations describe the local evolution of the large-scale structures 
in the solution. Notably, in the case of (11. lap , these equations are of hyperbolic type. 

The solution of ( II. ip by 1ST is intimately connected with the spectrum of the Schrodinger 
operator, 

dx^ 

and the first step in Lax & Levermore's analysis was to replace the true spectrum with WKB 
approximations. This replacement creates a sequence of refiectionless potentials i/q'' (x) which 
converge to yo{x) in L^(]R) as e J, 0. Using the 1ST, Lax & Levermore were then able to write 
down essentially explicit representations of the family of solutions of (11. lap corresponding 
to this sequence of modified data, and they were able to analyze and describe the limiting 
structure of this family. 

These ideas and their extensions have also been used to address other problems including 
the semiclassical limit of the defocusing nonlinear Schrodinger equation [23], the continuum 
limit of the Toda lattice [S], and a continuum limit of a discrete nonlinear Schrodinger 
equation [ID]. In all of these analyses, a step corresponding to the modification of the initial 
data, as described above for (II. ip . has been the starting point of the analysis. In each of 
these cases, the Whitham equations are hyperbolic, hence locally well posed. Thus, in view 
of the L^- convergence of the modified data to the true data, it seems reasonable to expect 
that the convergence holds for t > as well. Here, however, we address a case in which the 
Whitham equations are elliptic, and such an expectation seems much more dubious. Our 
aim here is to better understand the effect of modifying the data in such a case. 

1.2. Focusing Nonlinear Schrodinger Equation and Semiclassical Limit. We con- 
sider the initial-value problem for the semiclassically-scaled focusing nonlinear Schrodinger 
(NLS) equation: 

iedtu + —d^u + \u\'^u = , (1.3a) 
u{x,0) = uo{x) . (1.3b) 

Equation (11.3aP is a universal model equation that arises in models of diverse physical sce- 
narios; it describes the envelope dynamics of a monochromatic wave in a weakly dispersive 
nonlinear medium in which diffusive effects are negligible [H ESI SI] • For example, it is a 
simple model for the propagation of light in optical fibers [S|. In (ll.3ap . < e ^ 1 is a 
constant parameter which measures the ratio of dispersion to nonlinearity. 

Our interest is the semiclassical or zero- dispersion limit of (ll.3p . That is, we suppose that 
the initial data uq is fixed, and we solve (II. 3p for each small e > 0. We describe below our 
assumptions on uq which guarantee the existence of a unique global solution to (II. 3p so that, 
in principle at least, this first step is possible. Then, given the resulting family (indexed by 
e) of solutions, 

u{x, t) = M*-^^(x, t) , 

the goal is to describe the asymptotic behavior of these solutions in the limit e J, 0. The first 
breakthrough for this problem was due to Kamvissis, McLaughlin, & Miller [2S] for initial 
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data of the form 

uo{x) = Ao{x) , (1.4) 

where : M — )■ (0, A] is even, bell-shaped, and real analytic. More precisely, Aq is assumed 
to 

(i) decay rapidly at ±00; 

(ii) be an even function, i.e., Ao{x) = Aq{—x) for all x G M; 

(iii) have a single genuine maximum at x = 0, i.e., Aq{0) = A, A'q{0) = 0, Aq{0) < 0; and 

(iv) be real- analytic. 

Henceforth, we adopt these assumptions. 

We remark that for fixed e > 0, well-posedness for the Cauchy problem (with e-independent 
data, as in fll.3b[) ) is well known. For example, we note that Ginibre & Velo [22] have shown 
that if uo e H'^iR) fl L°^(M), then fOaj) has a unique global solution u{t) in ^(M; H\R) n 
L°°(M)); the solution depends continuously on the data. Moreover, in the case (as we consider 
here) that mq € ^(M) — the Schwartz space of rapidly decaying functions, it is known that 
u{-,t) is also in ^(R) for each t [L9\. The issue is that this well-posedness is not uniform in 

em- 

1.3. Inverse scattering transform. As in the other problems to which Lax-Levermore 
theory has been applied, equation (11. Sal) is integrable, i.e., there is an associated Lax pair. 
Thus, we can solve the initial-value problem (II. 3p by 1ST; [A] contains an outline of the 
process. Indeed, it could be argued that the integrability of (ll.3ap is the only feature that 
makes the task of obtaining (postbreak) asymptotics even appear tractable. The equivalent 
problem for nonintegrable variants of (ll.Sap appears to be widely open |10j. 

The first step, then, in solving (II. 3p is an analysis of the nonself adjoint Zakharov-Shabat 
eigenvalue problem (one half of the Lax pair for (ll.3al) ): 



dx 



Wi{x; A) 
^2(3;; A) 



— iA Aq{x) 
-Aq{x) iA 



Wi{x; A) 
W2{x; A) 



1.5) 



In (II. Sp . Wi and W2 are auxiliary functions and A G C is a spectral parameter. For each e > 
and for Aq as described above, it is known (see |28j) that the discrete spectrum of (II. 5p is 
confined to the imaginary axis. Beyond this, a formal WKB method applied to (II. 5p suggests 
that the reflection coefficient is small beyond all orders and the imaginary eigenvalues are 
given by a quantization condition of Bohr-Sommerfeld type. Since precise information about 
the true scattering data (discrete eigenvalues of (11.51) and the reflection coefficient) is not 
known, a natural way forward — following Lax & Levermore — is to use the (formal) WKB 
scattering data in its place. For each small e > this procedure, neglecting the reflection 
coefficient and using the WKB eigenvalues, amounts to replacing the true initial data Aq 
with some other initial condition u'^q-' which depends on e and for which the WKB spectral 
data is the true spectral data. Because we neglect reflection, each solution of (ll.3ap with 
initial data Uq^ is an A^-soliton with ~ e~^. The collection of all these exact A^-soliton 
solutions of (ll.3ap (with N 00 and e J, 0) is called the semiclassical soliton ensemble 
(SSE) associated with Aq. 

The analysis of Kamvissis et al. is almost wholly focused on the inverse scattering step 
for SSEs. Now, there are at least two special cases for the which the spectral data is known 
exactly. For Aq{x) = Asech(x), Satsuma & Yajima [3S], after a clever transformation, have 
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shown how to write down exphcit formulae for the eigenvalues, proportionality constants, and 
reflection coefficient coming from (11. 5p . More recently, Tovbis & Venakides [13] introduced 
a special family of initial data (with a complex phase) for which the forward scattering 
problem can also be treated exactly. This family of data forms the foundation of the related 
work on the semiclassical limit by Tovbis, Venakides, & Zhou [H]. It is clearly of interest 
to obtain results for a much wider class of initial data, and a natural first step is to look at 
the bell-shaped data that generate SSEs. 

Strictly speaking, the analysis of Kamvissis et al. [26] describes the asymptotic behavior 
of such SSEs for t ^ 0. At t = there is the complementary result of Miller: 

Theorem 1 (Miller [3^). In the situation described above, there is a sequence {ejy)'^^-^ such 
that 

lim ejv = 0, (1.6) 

and such that for each x ^ there exists a such that 

\ut\x) - Ao{x)\ < K,e]i'-\ iV = 1,2,3, .. . (1.7) 

for all u > 0. 

As noted above, a fundamental issue that remains unresolved is to connect the asymptotics 
of the SSE to those of the true solution of the initial- value problem (II. 3p . That is. Theorem [1] 
shows that that SSE and Aq are asymptotically pointwise close at t = 0. However, (ll.3ap 
has modulational instabilities whose exponential growth rates become arbitrarily large in the 
semiclassical limit — the Whitham equations are elliptic. Thus, it is not possible to conclude 
from Theorem [T] that any member of the SSE and the corresponding true solution are close 
for any t > 0. To attack this difficulty, one could try to rigorously estimate the deviation 
of the WKB spectral data corresponding to Aq from Aq's true spectral data. With such 
eigenvalue-by-eigenvalue control in hand, one could then try to incorporate this information 
into the asymptotic analysis. This is the ongoing work of ^ . Our complementary goal in this 
paper is to better understand, by numerical experiment, the relationship between the SSE 
and Ao aX t = and between the SSE and the true solution u{x,t) for t > 0. In particular, 
our experiments support the following conjecture. 

Conjecture. For small times, despite the presence of modulational instability, the particular 
asymptotically small modification of the initial data used by Kamvissis et al. to generate a 
semiclassical soliton ensemble for bell-shaped data does not affect the limiting behavior. 

In the remainder of this paper, we describe the process we used to generate numerical 
evidence that supports this conjecture. Remarkably, our experiments show that the 0(e) 
rate of convergence of the modified data to the true data is propagated to positive times, 
including times after wave breaking. 

1.4. Plan. To aid the reader, we now outline the contents of the remainder of this paper. 
In Section |2] we describe the machinery from the theory of integrable systems necessary to 
obtain an A^-soliton solution of (ll.3ap with initial data of the form (II. 4p . That is, assuming 
that refiectionless scattering data for Zakharov-Shabat problem (II. 5p is known, we recall that 
the solution of equation (ll.3ap can be obtained by solving an appropriate Riemann-Hilbert 
problem. In this section, for initial data Aq{x) = exp(— x^), we also describe the computation 
of the WKB eigenvalues for (II. 5p . With high-precision approximations of these eigenvalues in 
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hand, we are able to use known techniques |37[ 133] to generate members of the corresponding 
SSE. Sections [3] and m are devoted the development and testing of the numerical methods we 
use for our eventual comparison between members of the Gaussian SSE and the (numerically 
computed) true evolution of (11.31) . We use the members of the Satsuma-Yajima ensemble, 
known exact A^-soliton solutions, to validate our numerical methods in the range of e and t 
that we consider here. Using the methods of Section [3], we report on the principal experiment 
of the paper in Sectional this is the aforementioned comparison of true evolution for various 
values of e with that of the corresponding members of the Gaussian SSE. Finally, Section [6] 
contains a discussion of our results. For example, we contrast our results with some examples 
in the hterature which suggest that when e is small, equation (ll.3ap is extremely sensitive to 
(rough) perturbations of the initial data. \M contains an outline of the features of the inverse 
scattering transform for (ll.3al) that are used in this paper. 



2. RiEMANN-HlLBERT AND WKB 



2.1. Riemann— Hilbert formulation. We take as our starting point the fact that every 
A^-soliton solution of the focusing NLS equation can be characterized as the solution of a 
meromorphic Riemann-Hilbert Problem (RHP) with no jumps. The solution of the RHP is 
a matrix- valued rational function of A G C; the solution depends on a set of discrete data — a 
collection of complex numbers in the upper-half plane 

{^N,0, ^N,l, ■ ■ ■ , ^N,N-l} , (2.1) 

nonzero constants 

{7n,o,7n,i,---,1n,n-i}, (2.2) 
and a choice of J = ±1. One seeks to solve the following problem. 

Riemann-Hilbert Problem 1. Find a 2 x 2 matrix-valued function m{X;x,t) with the 
following properties. 

(1) m(A;a:,t) — )■ I as A — )■ cxo. 

(2) vn{\;x,t) is a rational function of A with poles confined to the values X^^k and A^;.. 
At the singularities, 

w [ o" 
resA=A^rfc m(A) = lim m(A)o-i " 



A— !>A» 



a 



i-j 

2 



resA=A^^, m(A) = lim m{X)al 



1 - J 
2 



CAr,fc(x,t) 





1- J 

2 



(2.3) 
(2.4) 



Here, 



CN.k{x,t) := ( — 
\7fc 



N-l 

Y\ {^N,k — A 



N,n) 



n=0 



N-l 



exp 



2iJ{XN,kX + XnM' 



(2.5) 



and o"! is the Pauli matrix 



1 

1 
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Finally, once the solution of RHP [T] is found, one recovers an A^-soliton solution via the 
formula 

u{x, t) = 2i lim Ami2(A; x, t) . 

A— >oo 

As pointed out by Kamvissis et al. the RHP above can be recognized as a classical 
Pade multipoint interpolation problem. They used RHP [1] as the starting point of their 
analysis; as a first step they exchanged the meromorphic problem above for a sectionally 
holomorphic one. Then, as the result of a substantial amount of work, they were able to 
transform the sectionally holomorphic RHP to one which is amenable to the steepest-descent 
techniques of Deift & Zhou [16]. The results of this elaborate analysis are detailed asymptotic 
formulae for the small-e behavior. 

We proceed in a different fashion. After making a partial fractions ansatz, it is possible to 
reduce the solution of RHP [T] to the solution of an x linear system; see [23 E3] or lAl 
Then, given eigenvalues {Xnj}, constants {■Jnj}, and a pair (x, t) (these appear in the hnear 
system as parameters), we may recover the A^-soliton solution of fll.3al) at {x,t). Thus, to 
construct members of a SSE associated to initial data Aq, we need to compute the WKB 
eigenvalues of the Zakharov-Shabat problem with Aq appearing as potential. With these 
in hand, we may then turn to solving the poorly conditioned linear system that is born of 
RHP [U We describe the calculation of the WKB eigenvalues in the next section. 

2.2. The WKB Formulae. 

2.2.1. General Case: Bell-shaped Data. We begin by recalling the formulae for the WKB 
eigenvalues of (11. 5p : for more details see pS]. The basic object of interest is the density 
function 

p\v) ■■= - / , = / ^/Mxf+^dx, (2.6) 

defined for t] G (0,1^4) where x±{ri) are the two real turning points; see Figured! From p° 



-ir] 





A y = Ao(x) 











Figure 1. The turning points x±{r]). 



we obtain the function 

e\\):=-7i (2.7) 
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which measures the number of WKB eigenvalues on the imaginary axis between A and iA. 



We then define for = 1, 2, 3, . . . 



-| piA 1 POO 



Finally, the WKB eigenvalues A^v^fc are defined (there are of them for e^) by the formula 

riA 



p°{r]) dr] = eN (k + 



1 



k = 



TC 



(2.8) 
(2.9) 



Using the above formulae we may "simplify" the left-hand side: 



iA 



x+(r?) 



^/A^)(xy~+~rf dx 



dr] 



2 r^+{^N,k) 



IT 



Ao{xy + ~X%i,dx 



Therefore, writing X^^k = i^iv,/c for t^^k £ (0, A) C M, we desire to solve the equation 
Ao{xr-t%^,dx = ^(k + -\ , fc = 0,1,2,... iV-1. 



(2.10) 



In this case, the auxiliary scattering data (proportionality constants) are given by 

lN,k = {-l)'^'. (2.11) 

2.2.2. The Gaussian SSE. For our numerical experiments, we restrict ourselves to the Gauss- 
ian SSE. That is, from now on, we consider the problem (11. 3p with fixed initial data given 
by 



Then, from fl2.8 



Uo{x) = Ao{x) = e" 
1 



7lN 



e dx 



and formula f l2.10p becomes 

'•a;+{itjv_fe) 



/ 

Jo 



-2x2 



— f2 dr 



"n 
2N 



k + -], k = 0,l,2,...N-l, 



where x± are given by 

x±{it) = ±V-\nt. 
Equation fl2.10p thus reduces in this case to 



-t%^dx = ^{k + - 



2N 



k = 0,l,2,...N -1. 



(2.12) 
(2.13) 

(2.14) 
(2.15) 

(2.16) 



It is the solutions tjy^k, k = 0, . . . , N — 1 of f l2.16p together with the 7Ar,A:'s which will generate 
the exact A^-soliton solution of f ll.3ap . The collection of these solutions for G N is the 
Gaussian SSE. 
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Our first task is to solve f l2.16p to very high precision. With the numerically computed 
WKB spectral data in hand, we then use the numerical linear algebra routines of [37] and 
[33] to reconstruct via inverse scattering various members of the SSE at t = (and later 
times too). High precision knowledge of the spectral data is necessary due to fact that the 
solution is obtained by solving a poorly conditioned linear system [37] . We will then compare 
the numerical reconstructions of members of the Gaussian SSE at t = to the true initial 
data Aq = and with approximations to the evolution at later times. 

We now make a few comments about the solution of (12.161) . We performed these calcu- 
lations with 250-digit precision in Maple. However, our initial attempts to solve equation 
(I2.16P directly were unsuccessful, and we found it necessary to transform the problem to 
avoid difficulties with the root finder. In particular, if we define 




(2.17) 



K{x,t) 



then equation (I2.16P can be viewed as solving the single equation 

F{t) = constant . 

Two views of the graph of the function K appearing in the definition of F in (12.170 are shown 
in Figure [2l and we attributed the failure of the Newton solver to the square-root vanishing 
of K and its infiuence on F'{t) along the curve (a;+(it),t). To overcome the problem and 




Figure 2. Two pictures of the graph of K. Note that K vanishes along the 
curve defined by {x^{it),t) = (V— Int, t). 

eliminate the difficulty, we define 

e-^^' =ecos\i^w, (2.18) 

and we write F{t) as 

F{t) = / Ve2^'' - t2 dx + / Ve2^' - t2 dx (2.19) 

=:Fi(t) + Fn(t). 



Then, changing the integral in Fu via f l2.18p to an integral with respect to w, we obtain 

Fn{t)= V e2^' - 12 dx = - / dw. (2.20) 

■^iV^ ^-^0 J-iln(t2cosh2w;) 



It follows that fl2.16p may be rewritten as 

/•Iv/^i^ ^— ^— — 1 ^cosh-^ft-y) t^^^ sinh ^; tanh w 
•^0 ^ 2io ./-iln(t^,cosh2^/; 



= ^(fc + M . (2.21) 
2A^ V 2y ^ ' 

This is the equation we solve to high precision. We verified the 250-digit accuracy of the 
solutions of fl2.2ip using both Mathematica and Maple routines. 

3. Numerical Methods 

We introduce an implicit finite difference scheme to solve the initial value problem (11. 3p 
in this section. In Section H] we use an exact solution to illustrate the order of accuracy of 
the proposed method. We show that the temporal grid sizes used for the proposed method 
scale linearly with the spatial mesh refinement. We then use particular A^-soliton solutions, 
members of the S at suma- Yakima ensemble [SB] , which we obtain by the 1ST calculation, to 
validate the proposed method for small-e calculation. At the same time, we compare the 
proposed method with the well-known spectral split-step method and show that the proposed 
method is a suitable method for solving the focusing NLS in the semiclassical regime. We 
also investigate a filtering process that removes Fourier modes whose amplitudes are smaller 
than a given threshold for our calculations with small e. Finally, in Section |5l we compare 
numerical solutions of the proposed method with that of the Gaussian SSE for the focusing 
NLS. 

The focusing NLS equation (ll.3a|] we consider here can be easily rescaled into the standard 
cubic NLS 

idt*ij + dl^tp + 2\ij\^ij = , (3.1) 
with u(t,x) = y/eip{2t*, y/ex*). Equation (13. ip is completely integrable in the sense of 1ST 
and has a canonical Hamiltonian form. A spatial finite-difference semi-discretization of 
equation (13. ip proposed by Ablowitz and Ladik [3], 

.dV^ + (ipm+l - 2ll)m + V^m-l) + IV^rn^ (V'm+l + V'm-l) = , (3.2) 

at 

is also completely integrable and posseses a Hamiltonian structure; here, h is the spatial 
grid size. We refer to the above discretization as the AL-lattice. Fornberg [20] has shown 
that with accurate (exact) time integration, the AL-lattice is very suitable for numerical 
work, since it produces few numerical artifacts for unstable analytical solutions in a periodic 
domain. More discussion on numerical homoclinic instability for the standard NLS can be 
seen, for example, in the papers by Ablowitz et al. [H [2]. Nevertheless, choosing a proper 
numerical time integrator for the AL-lattice is by no means a trivial task. 

Schober et al. [391 [23] indicate that the Hamiltonian system of the AL-lattice carries on its 
phase space a noncanonical symplectic structure for which standard symplectic integrators, 
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such as symplectic implicit Runge-Kutta methods, are not immediately applicable. Several 
approaches are provided by Schober et al. [391 123] to remedy the situation. While symplectic 
algorithms have the advantage of preservation of the global and local conservation laws for a 
long period of time, our aim is in the direction of developing an efficient and stable algorithm 
that is suitable for fine-grid calculations in order to accurately capture the behavior of fll.3al) 
in the small-e regime. 

Taking advantage of the simple and clean form of the AL-lattice, we propose to apply the 
implicit midpoint time integrator to the AL-lattice directly. The implicit midpoint method 
is the lowest order member of the Gauss-Legendre family of implicit Runge-Kutta methods 
which are symplectic schemes for canonical Hamiltonian systems [23]. Our numerical ex- 
periments show that the combination of the midpoint time integrator and the AL-lattice is 
advantageous for solving the semiclassical focusing NLS equation. The advantages include 
(1) the ratio of temporal grid size and the spatial grid size used for the method can be kept 
constant when refining the mesh; (2) with a good initial guess, the simple fixed-point iteration 
process converges relatively fast (< 10 iterations for the convergence tolerance 7 < 10~^^); 
(3) unlike the standard spectral split-step method, the proposed method is less sensitive to 
what spatial and temporal grid sizes to use in the simulations to avoid numerical artifacts 
caused by numerical roundoff error for small e. 

Finally, we remark that many numerical methods for the focusing NLS in the semiclassical 
regime have been discussed in the literature [El El [TJ [11] , but none directly compared with 
the 1ST calculation. 

3.1. Implicit Finite Difference Algorithm. The proposed finite difference scheme for 
the initial value problem (11.31) is as follows. 

Step 1. Based on the AL-lattice, the spatial discretization of equation fll.3ap is 
dum 1 

ie— TT^ + 777-^ (Wm+l - + Um-l) + -\Um\^ {Um+1 + Mm-l) = , (3.3) 

at 2Ax"^ 2 
where Ax is the spatial grid size, and Um approximates the solution at the m*'* grid point. 

Step 2. Applying the midpoint time integrator to the above ordinary differential equations 
(ODE) yields 

^n+l _ , '^^^t ("u^+^Z^ - 2m"+1/2 , y"•+l/2^ , i^|^n+l/2|2 (n+1/2 n+l/2\ .3 

where At is the time step size, and u^^^"^ is defined as 

«r = \K + . (3.5) 



Step 3. For n = 1 . . . A^, we solve the nonlinear equations f l3.4p by using the simple fixed- 
point-iteration (FPI) procedure, in which the {k + I)*'* iteration is written as 

n+Uk+l) n , ig^^ / n+l/2,(fc) _ ^ „+l/2,(fc) , n+l/2,(fc)\ 

-L l^|„n+l/2,(fc)|2 / n+l/2,(fc) n+l/2,(fc) 
~^ 2e y^m+l "T "m-l 
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(3.6) 



where u^'^^'^'^^^ is defined as 

<+V2,w = + . (3.7) 

The initial guess Mm^^'*^°'' for the FPI procedure within each time step is the solution of the 
Crank-Nicolson-type scheme for the AL-lattice: 



where u^^^'^ is defined in (13 .7^ . Equation (13 .Sp results in a tridiagonal system for u^^, which 
is solved by the Thomas Algorithm [27] . The convergence tolerance for the FPI procedure is 

||^n+l,(m)_^n+l,(fc)||^<^^ (3.9) 

where 7 < 10^^^ for the numerical experiments throughout this paper. Here || ■ ||oo is the 
infinity-norm defined by 

IImIIoo = max \um\ ■ (3.10) 

" m=l,...M 

When the convergence tolerance is achieved, we set = (*-■+!) ^ and move onto the 
next time step. 

3.2. Spectral Split-Step Method. Splitting schemes are very appealing for solving the 
focusing NLS equation in periodic domains. Within one At, a splitting method advances 
the NLS equation (ll.3al) by solving the following two equations alternately. 

(A) Nonlinear part (solve exactly in physical space) 

2i, ,2 , , 

Ut = —\u\ u . (3-11) 
e 

(B) Linear part (solve exactly in Fourier space): 

ut = \eu^x- (3.12) 

Yoshida [17] introduced a systematic method to construct arbitrary even-order time accurate 
splitting schemes. For example, to obtain second-order accuracy in time, we solve the two 
equations sequentially, like (A) — )■ (B) — )■ (A), by using the time increments {^, ^} 
in each step, respectively. Alternatively, one can also solve the sequence (B) — )■ (A) — > 
(B) to obtain the same order of accuracy, although this sequence is more time consuming, 
since one has to compute the (inverse) Fast Fourier Transform twice. We remark that the 
second-order accurate method constructed by using the Yoshida's scheme is essentially the 
Strang splitting method [Tj. The sequences and time increments for fourth and sixth-order 
methods are listed, for example, in the paper by Fornberg & DriscoU |21j . 

In Step (A), we solve the ODE (13.111) exactly. Bao et al. [7| did a simple calculation to 
show that |up in equation (13. lip is invariant within each time increment, 

4 4 

dt\u\^ = 2Re{utu) = -Re{i\u\\u) = -Re(i|M|^) =0. (3.13) 

e e 



Equation fl3.13p implies that the ODE (13. lip is linear and separable. 
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We use the second-order scheme to illustrate the split step algorithm as follows. Suppose 
that in a periodic domain, equation (11 .Sap is solved on a mesh for Um, where m = 1 . . . M + 1 
and um+1 = Ui, within a time step At = t"-~^^ — t". 

Step 1: We solve (A) in the physical space with the initial time and the final time 
t" + \At. The solution at the m*'* grid point is 

+ ^At) = M„(r)e'i'l"™(*")l' . (3.14) 

Step 2: We solve (B) in the Fourier space. The Fourier transform of equation (I3.12p is 

Ut = —ik'^eu , (3.15) 

where k is the wavenumber, and u is the Fourier transform of u. The Fourier transform pair 
for u are defined as 

M 

E(m-l){k~l) 

m=l 



M 

i 

Ur 



(3.16) 

k=l 

where 

a;^, = e(-2-)/^^ (3.17) 

Solving the ODE (I3.15P with the initial time + and the final time t" + |At yields the 
fc*'* Fourier mode of u, 

ur,{e + ^At) = uuit^ + iAt)e-5i'='^^* . (3.18) 

Step 3: After taking the inverse Fourier transform of u in f l3.18p . we solve (A) again in the 
physical space with the initial time + |At and the final time f^^^. The solution at the 
m*^ grid point is 

«^(r+i) = „^(r + ^At)e'i'l"™(*"+t^*)l' . (3.19) 

4. Numerical Experiments 

4.1. Exact Solution. It is easy to check that an exact solution associated with the focusing 
NLS equation (ll.3ap is 

u{x, t) = sech ^-^i-+'it) (41) 

The structure of the solution (14.11) is simple; it features a bell-shaped envelope propagating at 
a constant rate. We use this exact solution to validate our numerical implementation and to 
test the order of accuracy of the numerical methods. Tabled] is the mesh refinement study of 
the proposed implicit finite difference method and the split-step method. In the calculations, 
the length of the periodic domain is [—16, 16], the small parameter is e = 0.5, and the final 
time is t = 0.5. The temporal grid size for the finite difference method is At/ Ax = 1/80, 
and the convergence tolerance is 7 = 10~^^. If we define the physical quantity 

p{x,t) = \u{x,t)\\ (4.2) 
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the error of p between the numerical calculation and the exact solution measured by the 
2-norm 



1 



,2 

m ' 



(4.3) 



m=l 



is shown in Table [H The table indicates that both methods are second-order accurate. We 
also note that for both methods, At/ Ax can be kept roughly constant throughout the mesh 
refinement study to maintain the desired accuracy. 

Table 1. Mesh refinement study 





1/32 


1/64 


1/128 


1/256 


1/512 


1/1024 


IIPFD — Pcxact 2 


5.0836e-04 


1.2709e-04 


3.1772e-05 


7.9429e-6 


1.9857e-6 


4.9643e-07 


At /Ax 


1/80 


1/80 


1/80 


1/80 


1/80 


1/80 


order 




2 


2 


2 


2 


2 


IIPSS — Pexactlh 


1.1046e-06 


2.7624e-07 


6.8953e-08 


1.7372e-8 


4.4907e-9 


1.1117e-09 


At /Ax 


1/10 


1/10 


1/10 


1/10 


1/10 


1/13 


order 




2 


2 


1.99 


1.95 


2.01 



4.2. A^-soliton. Lyng & Miller |33] introduced accurate numerical reconstructions of the 
A^-soliton by the 1ST; these refined earlier calculations of Miller & Kamvissis [37]. In this 
case, the iV-soliton is the solution of the initial-value problem 

iV'* + ^V'xx + |^|V = 0, ^^^^ 
ipi^x, 0) = N sech(a;) . 

If the amplitude and the time variable of are scaled by a parameter e for a new variable, 
u{x, t) = eip^x, t/e), the initial value problem (14.41) is equivalent to the focusing NLS equation 
for u: 

ieut + ^e^Uxx + \u\'^u = , , , ^, 

2 ' ' (4.5) 

m(x, 0) = y4 sech(a;) , 

where e = A/N . It is well-known that the A^-soliton breaks its focusing state into the 
oscillatory state at the critical time tc = {2A)~^ We consider the case of initial amplitude 
A = 2 (thus tc = 0.25) and compute the initial value problem fl4.5p for various N by 
using the 1ST. We choose the time slice at t = 0.3 (> tc) with 4096 points in the interval 
< X < 1 on the {t, x) plane. All calculations are done using Mathematica with 250-digit 
precision. The solutions obtained by the 1ST calculation are, up to the numerical precision 
of the implementation of the 1ST, exact solutions of the A^-soliton problem (14. 5p . We use 
these solutions to assess the performance of the finite difference and the spectral split-step 
algorithms for small e. 

We now test these two numerical algorithms for the A^-soliton problem with the initial 
data u{x,0) = 2sech(x). The periodic computational domain is —16 < x < 16. We note 
that results obtained by using the initial data u{x, 0) = 2 sech(x) (e = 2/N) at the final time 
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t = tf are equivalent to those that computed with initial data u{x,0) = sech(a;) (e = 
at the final time t = 2tf, modulo a factor two in amplitude. 

Example 1: u{x,0) = 2sech(x), = 40, e = 0.05, final time t = 0.3. 

Figure [3] shows simulations for = 40 at the final time t = 0.3. Figure El^a) plots the 
conserved quantity p computed by using the finite difference method. Figure El^b) is the 
counterpart of Figure [3]|^a), using the spectral split-step method. Both (a) and (b) are 
plotted against the 1ST solution. We observe that numerical solutions of both methods are 
visually indistinguishable from the quasi-exact solution. For the finite difference method, 
the 2-norm error against the 1ST solution is 9 x 10"'^, and for the spectral split-step method, 
the error is 7.5 x 10^'^. All 2-norm errors in this section are measured for —1 < a; < 1, unless 
specified otherwise. The finite difference method uses the spatial grid size Ax = 1/4096 
and the temporal grid size At/ Ax = 1/300. The convergence tolerance is 7 = 10~^^. The 
spectral split-step method uses Ax = 1/4096 and At/ Ax = 1/10. Similar to the example 
with the exact solution f l4.ip . this experiment indicates that the spectral split-step method 
is able to capture the right solution, and is more efficient than the finite difference method 
for e = 0.05. We note that similar output to Figure [3] can be obtained by using e = 0.025 
with the initial data u{0,x) = sech(x) at the final time t = 0.6. We choose A = 2 for a 
shorter breaking time. 



N-40, dx=1/4096. dt/dx=1/300 N-40, dx=1/4096, dt/dx=1/10 




Figure 3. A^ = 40 at t = 0.3. (a) A comparison of the finite difference 
method and the 1ST. The 2-norm error is 9 x 10~^. (b) A comparison of the 
spectral split-step method and the 1ST. The 2-norm error is 7.5 x 10~^. 



Example 2: u{x,0) = 2sech(x), A^ = 50, e = 0.04, final time t = 0.3. 

Figure m shows a refinement study of the finite difference method for the 2sech(x) initial 
data with A^ = 50. The figure plots the computed quantity p against the 1ST results. 
Figure m^a) uses Ax = 1/2048 and At/ Ax = 1/300, and Figure Sl^b) uses Ax = 1/4096 and 
At/ Ax = 1/300. We observe that the 2-norm errors for (a) and (b) are 0.1519 and 0.1300, 
respectively. The convergence tolerances are 7 = 10~^^ for Ax = 1/2048 and 7 = 10"^^ for 
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Ax = 1/4096. The study suggests that for the finite difference method, combined spatial 
and temporal refinement will reduce the error for e = 0.04. We now compare Figure 111(b) 
with the previous example, Figure |3]^a). We observe that while these two calculations use 
the same spatial and temporal discretization, the error in Figure ID^b) is one order larger 
than that in Figure [31(a). The only difference between these two simulations is the small 
parameter e, for which one is 0.05 and the other is 0.04. This seems to indicate that when 
the focusing NLS equation becomes just a little bit more singular (i.e. e decreases from 0.05 
to 0.04), the roundoff error sets in rather swiftly. The phenomenon of roundoff error becomes 
much more prominent for the spectral split-step method. 



N-50, dx=1/2048, dt/dx-1/300 N=50, dx-1/4096, dt/dx-1/300 




Figure 4. = 50 at t = 0.3. Refinement study of the finite difference 
method, (a) As = 1/2048, At/ Ax = 1/300, and the 2-norm error is 0.1519. 
(b) Ax = 1/4096, At/ Ax = 1/300, and the 2-norm error is 0.1300. 



Similar to Figure HJ Figure [51 is a refinement study of the spectral split-step method. In 
this study, we explore how temporal and spatial grid sizes used for the spectral split-step 
method affect the simulations. We also investigate a filtering process introduced by Krasny 
[29] for ill-posed initial value problems, such as the vortex sheet roll-up problem. Krasny's 
experience with vortex sheet simulations suggests that if the problem is ill-posed, fewer grid 
points should be used for the simulation. The reason for this is that more grid points in- 
troduce shorter wavelengths into the numerical solution, and once the short wavelengths are 
spuriously perturbed by roundoff error, the computation collapses quickly. Krasny proposed 
a filtering process, now known as the Krasny filter, that eliminated Fourier modes whose 
amplitudes were smaller than a threshold to restore smoothness of the roll-up. For the focus- 
ing NLS equation (II. 3p . it is known that the problem becomes ill-posed when e approaches 
to zero [11]. Hence, the Krasny filter is sometimes applied to simulations of the equation in 
the semiclassical regime [71 [TT] . In particular, Bao et al. [7j showed that the Krasny filter 
successfully restored symmetry for simulations that showed breaking of symmetry. 

In the recent papers of Bao et al. |7J and Jin et al. [25], the authors suggested that 
typically the temporal and spatial grids used for the focusing NLS equation should satisfy 

Ax = 0(e), At = o(e). (4.6) 

Taking this suggestion and the hint from our numerical experiment for = 40 in Figure [2]^b), 
we first use the same spatial and temporal grid sizes. Ax = 1/4096 and At = (l/10)Ax = 
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1/40960, as that used in Figure [3](b) for = 50. On the left-hand-side of Figure [5]^c), we 
show that using this set of mesh sizes does not produce a satisfactory result, compared with 
the 1ST calculation. We apply the Krasny filter to the same calculation so that if \uk\ < rj, 
where \uk\ is the amplitude of Uk and Uk is defined in equation (13.161) . we manually set Uk 
zero. Here t] is the threshold level. On the right-hand-side of Figure [S](c), we show that 
when the Krasny filter [i] = 10~^^) is applied at the end of each time step to the simulation, 
the filtered result has a better match with the 1ST calculation. The 2-norm error is reduced 
to 1.1177 from 2.6949, and symmetry has mostly been restored. Our numerical experiments 
show that if we use finer temporal grid sizes, with or without the filtering process, the results 
only get worse. If we, on the other hand, decrease the temporal grid size, by trial-and-error, 
we find that when At = 2x 10~^, without the filtering process, we obtain a reasonable match 
with the 1ST result, as shown on the left-hand-side of Figure [5](b). If we keep this temporal 
grid size and coarsen the spatial grid size to Ax = 1/2048, without the filtering process, 
we obtain an even better match with the 1ST calculation, as shown on the left- hand- side of 
Figure [S](a). On the right-hand side of Figure O^a) and (b), we show that the Krasny filter 
does not improve the results. On the contrary, more significant phase-shift is observed for 
calculations with the Krasny filter. We note that we obtain almost identical results for the 
threshold levels from 10^^^ to 10^^°. In the next numerical experiment, we further show 
that the Krasny filters not only do not improve the results, they produce results that do not 
match with the 1ST calculations for small e. 

Example 3: m(x,0) = 2sech(x), = 54, e = 1/27 ^ 0.037, final time t = 0.3. 

In this example, we show that for = 54 (e ~ 0.037), the proposed finite difference 
method captures the proper waveform of the solution, while the solution obtained by using 
the spectral split-step method is heavily influenced by numerical artifacts. We further show 
that the Krasny filter not only fails to reduce the numerical artifacts in both methods, but 
produces solutions that are drastically different from the 1ST calculation. Figure [6t^a) shows 
that without the Krasny filter, the finite difference method produces a solution that match 
the 1ST result closely. We note that the spatial grid size is Ax = 1/2028 in this simulation. 
Using either a finer grid such as Ax = 1/4096 or a coarser grid such as Ax = 1/1024 will 
degenerate the result, regardless the choice of temporal grid size, based on our numerical 
experiments. The convergence tolerance is 7 = 10~^^ for the simulation. Figure MJo) shows 
the same computation, except Krasny's filter [i] = 10~^^ and 10~^^) is applied at the end of 
each time step. The numerical results are identical for these two different thresholds, and 
they fail to match the 1ST calculation. 

Figure [7](a) shows that the result by using the spectral split-step method does not match 
the 1ST calculation for N = 54. The spatial grid size is Ax = 1/1024 in this simulation. 
Similar to the finite difference method, using either a finer grid such as Ax = 1/2048 or 
a coarser grid such as Ax = 1/512 will not improve the result. When the Krasny filter is 
applied to this simulation, the result becomes even worse, as shown in Figure [7](b) . Similar 
to the finite difference method, two thresholds, rj = 10~^^ and 10~^^, are applied to the 
computation, and the outcomes are identical. 
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N=50, dx-1/2048, dt-2e-4 



N-50, dx=1/2048, dt=2e-4, Ti=10"^ 




Figure 5. = 50 at t = 0.3. Refinement study of tlie spectral split-step 
method. Left column: without the Krasny filter. Right column: with the 
Krasny filter and the threshold is r/ = lO'^^. (a) Ax = 1/2048, At = 2 x lO""^. 
The 2-norm error is 0.4397 without the filter and 0.7561 with the filter, (b) 
Ax = 1/4096, At = 2 X 10~^. The 2-norm error is 0.5527 without the filter and 
0.7857 with the filter, (c) Ax = 1/4096, At = (1/10) Ax = 2.44140625 x 10~^ 
The 2-norm error is 2.6949 without the filter and 1.177 with the filter. 
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The results of this numerical experiment are hardly surprising. As pointed out by Jin et 
al. [25] , the Krasny filtering process actually violates the conservation of mass, which means 
that the cut small-amplitude Fourier modes could very much be part of the solution. 



N-54, dx=1/2048. dt/dx=1/200 N=54, dx-1/2048, dt/dx-1/200 




-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 -0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



(a) ' (b) 

Figure 6. N = 54, Ax = 1/2048, At/ Ax = 1/200, t = 0.3. Computation by 
using the finite difference method, (a) Without the Krasny filter, (b) With 
the Krasny filter. The thresholds are rj = 10^^^ and 10~^^, respectively. 



N=54, dx-1/1024, dt-2.5e-4 



N-54, dx=1/1024, dfc2.5e-4 
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-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



(b) 

Figure 7. = 54, Ax = 1/1024, At = 1200, t = 0.3. Computation by using 
the spectral split-step method, (a) Without the Krasny filter, (b) With the 
Krasny filter. The thresholds are rj = 10~^^ and 10~^^, respectively. 



5. Evolution: Aq versus u^^ 

In this section, we compare numerical solutions obtained by the evolutionary numerical 
methods discussed in Section [3] and by the 1ST for the Gaussian SSE described in Section [2l2l 
First, for the initial data given in equation fl2.12p . we compare the numerical solutions ob- 
tained by the finite difference method and the split-step method to benchmark the solutions 
for this initial data. Table |2] shows the 2-norm difference of p between these two methods 
at final times t = 0.1, 0.2, 0.3, 0.4, and 0.5, with number of solitons, = 5, 10, 15, and 

18 



20. We recall that the relation between the small parameter e and the number of solitons 
for this Gaussian initial data is given in equation f l2.13p . The table shows that at this 
range of e and final times, numerical solutions obtained by these two methods follow each 
other closely. The computational domain is — 10 < x < 10 and periodic. The grid resolution 
is Ax = 1/4096 for both methods, and the difference is measured for — 1 < x < 1. Figure 
[H] shows the case that has the largest difference in Table [21 for which = 20 and the final 
time t = 0.5. The two results are visually indistinguishable under the scale as shown in (a), 
and there is a moderate visual discrepancy after we magnify the graph, as shown in (b), for 
the centered peak. We remark that our numerical investigation for the Gaussian initial data 
below is limited to these ranges of and t. 

Table 2. 2-norm difference of p between the finite difference and the split- 
step methods for the Gaussian initial data. 





0.1 0.2 0.3 0.4 0.5 


5 


3.8360E-8 1.8054E-7 5.8820E-7 2.3012E-6 1.9354E-5 


10 


1.6094E-7 7.7978E-7 2.7336E-6 2.0142E-5 1.1714E-4 


15 


3.6547E-7 1.7772E-6 6.3417E-6 7.5820E-5 5.6675E-4 


20 


6.5206E-7 3.1748E-6 1.1401E-5 2.1776E-4 1.7275E-3 



Gaussian N=20, t=0.5, dx=1/4096 



Gaussian N-20, t-0.5, dx-1/4096 




-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 




(b) 

Figure 8. Gaussian data for = 20 (e 0.0282) at final time t = 0.5. The 
grid resolution is Ax = 1/4096. (a) Comparison between the finite difference 
method and the split-step method, (b) Magnification of the centered peak in 
fa). 



Figure |9] shows the reconstruction of the initial data associated with the Gaussian SSE, 
and its comparison with the true Gaussian initial data for A^ = 5 and 20, respectively. These 
reconstructed initial data can be seen as perturbation of the true Gaussian data. Our main 
interest is to understand how these perturbed data evolve with time in comparison to the 
evolution of the true data. For large A^ (small e), while the initial perturbation is small, how 
the perturbed data will evolve with time in comparison with the true evolution is not clear. 
In principle, the evolution of the perturbed data could depart from the true evolution in a 
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short time, if e is small enough. Figure [TO] shows that at time t = 0.5, the perturbed data 
become visually indistinguishable from the finite-differenced solution for both N = 5 and 
= 20. The computational domain is —10 < x < 10 and periodic, with the grid resolution 
Ax = 1/4096 for the finite difference method. 




(a) ' (b) 

Figure 9. Reconstruction of the Gaussian initial data by the Gaussian SSE. 
(a) N = 5,{h)N = 20. 




-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 -0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 



(a) ' (b) 

Figure 10. Comparison of a member of the Gaussian SSE (computed via 
1ST) with the corresponding finite-difference solution with Gaussian initial 
data at t = 0.5. (a) = 5, (b) = 20. 

We now systematically compute the 2-norm difference of p between the SSE solutions 
and those of the finite difference method (representing the true evolution of the initial data 
Mo = exp(— x^)) for the ranges of = 5 to = 20 and t = 0.0 to t = 0.5. Table |3] shows 
that the differences between the two solutions diminish with time for all A^ (except the last 
point, when A^ = 20, t = 0.5). Figure [TT] represents the same data; the markers show the 
2-norm differences versus A^ for times t = 0.0, 0.1, 0.2, 0.3, 0.4, and 0.5, respectively. The 
figure shows a remarkable consistency in the decay of the error as A^ increases; we discuss 
this further in the next section. We also note that the 2-norm difference between the finite 
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Table 3. 2-norm difference of p between the finite-difference solution (for t > 
0) with Gaussian initial data and the corresponding member of the Gaussian 
SSE. 



N 


0.1 0.2 0.3 0.4 0.5 




3 47SfiE-9 3 2Qfi5E-2 9 7883E-9 9 0fi79E-9 1 3330E-9 8 8719E-3 




9 Q90QE-9 9 7575E-9 9 3040E-9 1 fifiQ7E-9 1 0454E-9 8 7774E-3 


7 


9 5037E-9 9 3543E-9 1 Q490E-9 1 3797E-9 8 '^44QF,-3 5 Q334E-3 


8 

(J 


2 1881E-9 2 04Q8E-2 1 6701 E-9 1 1 51 7E-2 fi 7fi53E-3 5 34fi7E-3 

. X (J (J X XL/ ^ . W^<L/(JXL/ X . \J 1 W X XL/ ^ X . X 5j X 1 XL/ \J . 1 \J tJOXL/ O tj . 0^\J 1 X_/ O 


Q 


1 Q4Q0E-9 1 8900E-9 1 4fi77E-9 Q Q944F,-3 5 71 fi3E-3 4 4n58E-3 

X . U^UVJ Hi jLi X .(J jLi\J\J1U X .t:\J I 1 Hi jLi O . Z(T:t:XL/ O . 1 X VJOXL/ O t:. t:W tJOXL/ O 


1 n 


1 7fi98E-9 1 fi41QE-9 1 313fiE-9 8 7fi82E-3 5 0452E-3 3 7535E-3 

X. 1 \j LjKJ Hi X.\J^XiyX_/ ^ X.OXOxJXL/ (J. ( \J(J.ijXL/ O 'O.Xj^xJ^^Hi kj O. ( xJxj'OHi xj 


1 1 

X X 


1 fiOQ5E-2 1 4Q5QE-9 1 1 88QE-2 7 8fil 5E-3 4 5fil 5E-3 4 0957E-3 

X .\j\jZ)0 Hi Li ± .rtc/O Hi X . X (J O C/ Hi Li 1 . (J VJ X tj Hi kj t: . tj VJ X tj XL/ kj rt.VJZjO I Hi O 


12 


1 4772E-2 1 3fiq8E-2 10812E-2 7 0770E-3 4 138QE-3 3 fi043E-3 


13 


1.3611E-2 1.2591E-2 9.8614E-3 6.3696E-3 3.7160E-3 3.5731E-3 


14 


1.2605E-2 1.1631E-2 9.0362E-3 5.7488E-3 3.3310E-3 2.6767E-3 


15 


1.1751E-2 1.0818E-2 8.3428E-3 5.2359E-3 3.0135E-3 2.8144E-3 


16 


1.1025E-2 1.0131E-2 7.7671E-3 4.8292E-3 2.7846E-3 2.4003E-3 


17 


1.0389E-2 9.5319E-3 7.2736E-3 4.4972E-3 2.6298E-3 2.3538E-3 


18 


9.8091E-3 8.9850E-3 6.8248E-3 4.1976E-3 2.4770E-3 2.2044E-3 


19 


9.2696E-3 8.4754E-3 6.4019E-3 3.9052E-3 2.3439E-3 1.9980E-3 


20 


8.7755E-3 8.0077E-3 6.0095E-3 3.6229E-3 2.1483E-3 2.7102E-3 



difference and the split-step method is of the order 10~'^ for = 20 and t = 0.5, and this 
is also the difference between the finite difference method and the 1ST. Therefore, at this 
point it is difficult to determine whether the difference between the finite difference method 
and the 1ST truly represents the difference between the true solution and the IST-generated 
member of the SSE. Also, it has little meaning to do more comparison for larger t or N 
beyond this point. 

6. Discussion and Concluding Remarks 

We have compared a spectral split-step method and an implicit finite difference method 
for solving the focusing NLS equation in the semiclassical regime. In the special case that 
the initial data is Aq{x) = y4sech(x), the 1ST solution serves as an exact solution for the 
comparison. We find that the spectral split-step method is more efficient compared with the 
proposed implicit finite difference method. However, for small e, we find that the proposed 
implicit finite difference method is less sensitive to the choice of spatial and temporal grid 
sizes than the split-step method; poor choices lead to numerical artifacts caused by numerical 
roundoff error. We observe that to obtain simulations with the fewest numerical artifacts 
for the iV-soliton problem with large (e.g. A^ > 54 and A = 2), the use of spatial and 
temporal grid sizes for the numerical methods (both split-step method and the proposed 
finite difference method) should follow Krasny's suggestions ^29(1, i.e., use fewer grid points 
and larger time-step sizes, rather than the meshing strategy in equation fl4.6p . 

We also investigated a filtering process, known as the Krasny filter. We find that the 
process may help to restore symmetry of the solution, but the restored solution may not 
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Figure 11. The 2-norm differences versus for t = 0.0, 0.1, 0.2, 0.3, 0.4, 
and 0.5. The data points are the computed values of the 2-norm error E = 
IIpfd — Psse||2 from Table [31 The plotted curves are of the form E = C ■ A^" 
where the constants C, a are determined by a least squares fit to the data 
points; see Table HI The values for a in the legend show a 0(l/iV) = 0(e) 
rate of convergence even for the times t = 0.4 and t = 0.5 which are after the 
breaking time. 

represent the solution of the problem. Furthermore, when e is small, the filtering process 
could even destroy good numerical simulations of the problem. This is because for small e, 
the highly oscillatory analytical solution could be a superposition of those small-amplitude 
Fourier modes, it is not possible for the filter to distinguish the small-amplitude Fourier 
modes of the solution from those due to roundoff error. 

Finally, we used the two studied numerical methods to investigate the Gaussian SSE for 
the focusing NLS equation. Within the range of e and t for which we are confident about 
our numerical solutions, we find that for larger e, the perturbation of the initial data is 
quickly dissipated and the SSE solution becomes close to the finite differenced solution, 
a good approximation to the true solution. We see this as a reflection of the particularly 
special nature of the perturbations we consider here; after all, they are connected to the data 
through the WKB analysis of f ll.5p . By way of comparison, we recall one of the experiments of 
Bronski & Kutz [S] which featured a non-analytic perturbation of initial data. In particular, 
they considered a perturbation, see Figure [T2| of the initial data Uo{x) = sech(a;) by a small 
multiple of the tent function 



In this case, Bronski & Kutz found that the evolution in their simulations was extremely 
sensitive to even small amounts of non-analyticity. Similarly, Clarke & Miller [T3| also found 
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Figure 12. The non- analytic perturbation of Bronski &; Kutz [S]. 
Table 4. The values of C and a for each t. 



t 


0.0 


0.1 


0.2 


0.3 


0.4 


0.5 


a 


-0.9922 


-1.0196 


-1.1063 


-1.2548 


-1.3016 


-1.0580 


C 


0.1727 


0.1713 


0.1671 


0.1574 


0.1044 


0.0489 



that the problem is quite sensitive to non-analytic perturbations; they found wild behavior 
even for small perturbations of class C^(M). 

To reveal further structure in the data assembled in Table [31 we postulate a relationship 
for the 2-norm error E as function of N of the form 

E = IIpfd — PssElb = C ■ N"" , 

and we use the data compute values for C and a by least squares. The resulting values are 
shown in Table HI and the corresponding curves are plotted in Figure [TH The values for 
a in Table H] show that the modified initial data arising from the WKB analysis of (11.51) 
converges to the true data at a 0{1/N) rate. (Recalling, (12.131) . we see equivalently, that 
this rate of convergence is 0(e).) Remarkably, our experiments show this rate of convergence 
is preserved at later times suggesting the possibility of a limited kind of well-posedness in 
the semiclassical limit. Also, we recall that Clarke & Miller [13] showed how to extract an 
upper bound for the breaking time — roughly, the first time that |M(-,t)p begins to exhibit 
e-scale oscillations — by an analysis of the formula 

\t\ = ^—\ dy, (6.2) 

where E{m) is a complete elliptic integral of the second kind, and Aq is the initial data. 
Using their technique, we compute the first critical point of t(po) and find an upper bound 
on the breaking time for Aq = exp(— x^) to be 

tub = 0.377417. 

Thus, our experiments show that the 0{1/N) rate of convergence that we see for small times 
persists even past wave breaking. This we see as especially noteworthy. Taken together, 
these results are consistent with the conclusion that, despite modulational instability, the 
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asymptotically small modification of the initial data used by Kamvissis et al. [26] does 
not affect the semiclassical limit. Still, it is difficult, based on our experiments, to draw 
definitive conclusions about the true limiting behavior. For our approach to yield more 
definitive results, better numerical schemes are required to further investigate the problem 
for smaller values of e and larger times. This is currently under our investigation. 
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Appendix A. Background: Inverse- Scattering Transform 



In this appendix, for the benefit of readers who may not be familiar with the inverse- 
scattering transform, we outline its relevant aspects here. Notably, we also outline the origin 
of the linear system that we solve to generate members of the Gaussian SSE; further details 
on these calculations can be found in [371133]. The procedure for solving the nonlinear initial- 
value problem (11.31) via 1ST is analogous to the procedure of solving initial- value problems for 
linear partial differential equations by Fourier transform. Briefly, the initial data is mapped 
to the scattering data; the scattering data have a simple evolution in time; and the solution 
at later times is reconstructed from the time-evolved scattering data. This discussion is 
modeled after [5B|. For more details, see, e.g., [H [T^ H5]. 

The Lax pair for the focusing NLS equation fll.3ap consists of the following two linear 
equations 



and 



dw 



dw 

dx 



-iA2 + 



-iA 

-u* 



u 
iA 



w 



=: Uw 



Am + T^dxU 



u\ 



w 



Vw. 



(A.l) 



(A.2) 



[U,V]=0 



(A.3) 



Evidently, the zero-curvature condition 

dt dx 

is equivalent to the NLS equation. Note that the left-hand side of ( \A.3\i is independent of 
A and vanishes exactly when u solves fll.3ap . It is precisely the existence of the Lax pair 
that allows the construction of a large family of exact solutions. Effectively, we are able to 
replace the nonlinear problem flL3ap with the pair of linear problems (lA.ip , (\A.2\i . 



A.l. Scattering. The first step is a careful study of the problem (lA.ip for A G M with 
u = Uo{x). The analysis is facilitated by the fact that |Mo(a;)| decays rapidly as |x| — )■ oo, 
whence U tends to the constant matrix —iXa^ in the limit. It is precisely this observation 
that allows one to construct particular solutions, the Jost solutions, of the linear system 
(lA.ip normalized at ±oo. The Jost solution matrices J±(x;A) (the 2x2 matrices whose 
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S(A) 



AgM. (A.4) 



columns are the normalized Jost solutions) are both nonsingular fundamental matrices for 
the differential equation. That is, 

ox 

Of course, the 2x2 system can only have two linearly independent column solutions, and 
therefore there is a 2 x 2 matrix S, called the scattering matrix, such that J+(x;A) = 
S(A)J_(x; A). Further analysis reveals that the scattering matrix can be written in the form 

^a(A)* 6(A)*" 
-6(A) a(A) 

Here, a and h are complex-valued functions, and they form the basis of the transmission 
coefficient T(A) = l/a(A) and the reflection coefficient -R(A) = 6(A)/a(A). 

Now, it turns out that the function a has an analytic continuation to the upper-half of 
the complex plane. Indeed, a careful look at the definition of a — it is a determinant whose 
columns are formed from Jost solutions one decaying at each of the spatial infinities — shows 
that zeros of a in the upper half plane are L^(M) eigenvalues of flA.l|) . Associated with 
each such eigenvalue A^, there is a complex number 7^. which is the ratio of the two analytic 
solutions of flA.ip which make up the Wronskian a. The refiection coefficient R does not 
generally extend off of the real line. 

A. 2. Time evolution. Now, if u{x,t) solves fll.3ap with initial data Uq{x), then for each 
t > the entries in the coefficient matrix U will change. This means that the eigenvalues 
{Afc}, the associated proportionality constants {7^} and the reffection coefficient could be 
computed independently for each positive time. However, equation f lA.2p constrains the 



temporal evolution, and it is possible to write down explicit formulae which describe the 
time evolution. In particular, the Jost matrices, which we now write as J±(x,t;A), must 
satisfy the differential equation 

^ = iA^J.aa + VJ, . 

This, in turn, is enough to derive a differential equation in t for the scattering matrix S(A; t): 

^ = iA^[S,a3]. (A.5) 

Writing f lA.SP in components, we immediately discover that when u{x,t) satisfies fll.3ap . the 
function a(A; t) = a(A) is independent of t, and 

It is immediate that the eigenvalues {A^} are independent of t and that the reflection evolves 
simply as -R(A; t) = R{X; 0)e^''*'^*. Finally, a brief calculation shows that 7A;(t) = 7fc(0)e^^'^^*. 

A. 3. Inverse scattering. The solution u{x,t) can be recovered from the scattering data 
(eigenvalues, proportionality constants, reflection). One way to visualize this process is to 
combine the columns of the Jost matrices to obtain matrices which extend into the half 
planes SJA ^ 0. These matrices are meromorphic functions of A on the disjoint half planes 
with (generically) simple poles at the A„'s and their complex conjugates. The residues at 
these poles can be computed. Moreover, the boundary values of these matrices do not 
generally match on the real line; their mismatch can be quantifled in terms of the reflection 
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coefficient. Finally, the large-A asymptotics of these matrices are prescribed. Moreover, the 
solution u of the NLS equation is encoded in the large-A behavior of the second column of 
this matrix. The inverse scattering process amounts to turning this process on its head; 
the properties of the matrix enumerated above (meromorphicity on Q'A ^ 0, prescribed 
poles and residues, prescribed jump across M, prescribed large-A behavior) are in most cases 
sufficient to determine the matrix itself. A convenient way to organize this information is in 
a Riemann-Hilbert problem as follows. 

Riemann— Hilbert Problem 2. Find a 2 x 2 matrix-valued function m{X;x,t) with the 
following properties. 

(1) m{X;x,t) is an analytic function on 

C \ (M U {Ao, Ai, . . . , Aat-i, Aq, . . . , X^i}) ■ 

(2) m(A;a;,t) — )■ I as A — > oo. 

(3) m has simple poles at the points A^ and A^; the residues satisfy: 



resA=Afe m(A) = lim m(A) 

A— s>Afe 



resA=A* rn(A) = lim m(A) 





ek{x,t) 

-ek{x,t) 




Here, ek{x,t) is given explicitly in terms of the 7fc's. 
(4) The matrix m{X;x,t) takes continuous boundary values on 



(A.6) 
(A.7) 

For A G M we write 



m±(A; X, t) := lim m(A ± i6; x, t) , 
54,0 



and 



m+(A; X, t) = m+(A; x, t)v(A; x, t) {K.i 
with the jump matrix v given explicitly in terms of the reflection R. 
Finally, once the solution of RHP [2] is found, one recovers the solution via the formula 

u{x, t) = 2i lim Ami2(A; x, t) . 



A— >oo 



A. 4. A^-soliton solutions, linear system. A solution for which -R(A) = is completely 
characterized by the N eigenvalues in the upper half plane (and their proportionality coeffi- 
cients). In this case the jump matrix satisfies v = I so there is no mismatch between m± on 
M. This is the case we consider in this paper. 
If we make the partial-fractions ansatz 



m 



(A;x,t) = i+x;^^+i;^'^'^''^ 



fc=0 



fc=0 



A -At 



it's clear from RHP [T]that 
and 



resA=Afc m(A; x, t) = Afc(x, t) 



resA=A* m(A; x, t) = Bfc(x, t) . 
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(A.9) 

(A.IO) 
(A.ll) 



Moreover, we see that the matrices Ak must have zeros in the second column, 

Ak{x,t) 



akix,t) 
bk{x,t) 



while the matrices have zeros in the first column 

Ck{x,t) 



Bk{x,t) 



Then 



Similarly, 







o" 


_bk_ 


= efc ^ 


1 


Ck 


= el( 


'l" 


_dk_ 








dkix,t) 
1 



N-l 

E 

j=0 

E 

j=0 



(A. - A*) 



(A. 12) 



(A.13) 



(A.14) 



(A.15) 



(A^ - A,) 

Evidently, f lA.14p . (lA.lSP form a linear system for the coefficients in and Bfc. 
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